if !d.name eq 'PS' then begin
    device,xsize=18,ysize=15,yoffset=3
    !p.thick=1 & !x.thick=1 & !y.thick=1
end
!p.charsize=1.4
!x.margin=[10,7]
; equador
lev=grange(omega_0-35,omega_0+60,nlev)
yr=[min(lev),max(lev)]

plot,[0.62,0.621],[yr[0],yr[0]],xr=[0.62,1],yr=yr, xtitle='!8r/R!6',$
     ytitle='!7X!6/2!7p!6 (nHz)',ystyle=1,xstyle=1

ieq=where(th GE -2. AND th LE 2.)
oplot,r,total(omega[*,ieq],2)/float((size(ieq))[1]),thick=5,col=48
;
; 15 degrees
i15=where(th GE 12. AND th LE 17.)
oplot,r,total(omega[*,i15],2)/float((size(i15))[1]),thick=5,li=2,col=80
;
; 30 degrees
i30=where(th GE 28. AND th LE 33.)
oplot,r,total(omega[*,i30],2)/float((size(i30))[1]),thick=5,li=3,col=112
;
; 45 degrees
i45=where(th GE 43. AND th LE 48.)
oplot,r,total(omega[*,i45],2)/float((size(i45))[1]),thick=5,li=4,col=176
;
; 60 degrees
i60=where(th GE 57. AND th LE 62.)
oplot,r,total(omega[*,i60],2)/float((size(i60))[1]),thick=5,li=5,col=180
oplot,[0.96,0.96],yr,li=2,thick=1
oplot,[0.718,0.718],yr,li=2,thick=1

; annotations
deg='!Uo!N'
annstart_x = 0.73
len_line = 0.04
len_line2 = 0.01
annstart_y = 460
delstart_y = 7
offset_y = 1
textstart_x = annstart_x+len_line+len_line2

xs_line = [annstart_x,annstart_x+len_line]
ys_line = [annstart_y,annstart_y]
xs_symb = [annstart_x+len_line2,annstart_x+len_line2]

oplot,xs_line,[annstart_y,annstart_y],li=0,thick=4,col=48
xyouts,textstart_x,annstart_y-offset_y,'0'+ deg
annstart_y = annstart_y-delstart_y

oplot,xs_line,[annstart_y,annstart_y],li=2,thick=4,col=80
xyouts,textstart_x,annstart_y-offset_y,'15'+ deg
annstart_y = annstart_y-delstart_y

oplot,xs_line,[annstart_y,annstart_y],li=3,thick=4,col=112
xyouts,textstart_x,annstart_y-offset_y,'30'+ deg
annstart_y = annstart_y-delstart_y

oplot,xs_line,[annstart_y,annstart_y],li=4,thick=4,col=176
xyouts,textstart_x,annstart_y-offset_y,'45'+ deg
annstart_y = annstart_y-delstart_y

oplot,xs_line,[annstart_y,annstart_y],li=5,thick=4,col=0
xyouts,textstart_x,annstart_y-offset_y,'60'+ deg
annstart_y = annstart_y-delstart_y

end

